***DRAW TABLE 4, City of 10,000 Cell

version 15.1

clear
use "worms_replication.dta"
set more off



*keep only covered dioceses*
drop if covered == 0


egen year_pop = group(year city_10_fixed)
collapse (mean) religious_bishop year city_10_fixed, by(year_pop)

rename religious_bishop avg_religious_p


*calculate and record the 40 year moving average share of religious bishops for dioceses with and without a city of 10,000 inhabitants before Concordats*
gen moving_avg_p = .
quietly forval i = 800/1299{
	sum avg_religious_p if year >= `i'-20 & year <= `i'+20 & city_10_fixed == 1
	replace moving_avg_p = r(mean) if year == `i' & city_10_fixed == 1
	sum avg_religious_p if year >= `i'-20 & year <= `i'+20 & city_10_fixed == 0 
	replace moving_avg_p = r(mean) if year == `i' & city_10_fixed == 0
	}

label var avg_religious_p "Share Religious Bishops"
label var moving_avg_p "Share Religious Bishops"
label var year "Year"

*For each year and each of the presence or absence of a city of 10,000 inhabitants plot both the share of religious bishops (avg_religious_p) and the moving average*
twoway (scatter avg_religious_p year if  city_10_fixed == 0  & year >=800 & year <= 1309, mcolor(black) lcolor(black))(line moving_avg_p year if city_10_fixed == 0  & year >=800 & year <= 1309, mcolor(black) lcolor(black))(scatter avg_religious_p year if  city_10_fixed == 1  & year >=800 & year <= 1309, mcolor(gs12) lcolor(gs12))(line moving_avg_p year if city_10_fixed == 1  & year >800 & year <= 1309, mcolor(gs12) lcolor(gs12) xline(1122, lpattern(dash) lcolor(black)) xline(1309, lpattern(dash) lcolor(black))), ///
   graphregion(color(white))  bgcolor(white) scheme(lean1) xlab(800 900 1000 1100 1200 1300)  ylab(, nogrid) legend(order(1 "poor diocese" 2 "poor diocese" 3 "wealthy diocese" 4 "wealthy diocese")) title("City over 10,000")

